#include <iostream>
#include <vector>
#include <cmath>
#include <complex>
#include <algorithm>
#include <memory>
#include <random>
#include <functional>
#include <fftw3.h>
#include <Eigen/Dense>

// 常用类型定义
using Complex = std::complex<float>;
using Pulse = std::vector<Complex>;
using PulseDopplerMatrix = std::vector<std::vector<Complex>>;
using Eigen::MatrixXcf;
using Eigen::VectorXcf;

// 雷达参数配置
struct RadarParameters
{
    float carrier_frequency;  // 载频 (Hz)
    float pulse_width;        // 脉冲宽度 (s)
    float pulse_rep_interval; // 脉冲重复间隔 (s)
    int num_pulses;           // 相参处理间隔内的脉冲数
    int num_samples;          // 每个脉冲的采样点数
    float sampling_rate;      // 采样率 (Hz)
    float bandwidth;          // 带宽 (Hz)
    int num_antennas;         // 天线阵元数量
    float antenna_spacing;    // 天线阵元间距 (m)
    float light_speed;        // 光速 (m/s)

    // 计算派生参数
    float getRangeResolution() const
    {
        return light_speed / (2.0f * bandwidth);
    }

    float getMaxUnambiguousRange() const
    {
        return light_speed * pulse_rep_interval / 2.0f;
    }

    float getMaxUnambiguousVelocity() const
    {
        return light_speed / (4.0f * carrier_frequency * pulse_rep_interval);
    }

    float getVelocityResolution() const
    {
        return light_speed / (2.0f * carrier_frequency * num_pulses * pulse_rep_interval);
    }
};

// 检测单元信息
struct DetectionCell
{
    int range_bin;
    int doppler_bin;
    int angle_bin;
    float power;
};

// 目标信息结构
struct TargetInfo
{
    float range;        // 距离 (m)
    float velocity;     // 速度 (m/s)
    float angle;        // 角度 (度)
    float power;        // 功率 (dB)
    float snr;          // 信噪比 (dB)
    int num_detections; // 凝聚的点迹数量

    TargetInfo(float r, float v, float a, float p, float s, int n)
        : range(r), velocity(v), angle(a), power(p), snr(s), num_detections(n) {}
};

// 信号处理工具类
class SignalProcessor
{
public:
    // 计算信号功率 (dB)
    static float computePowerDB(const Complex &sample)
    {
        return 10.0f * std::log10(std::norm(sample) + 1e-12f);
    }

    // 计算信号功率 (线性)
    static float computePowerLinear(const Complex &sample)
    {
        return std::norm(sample);
    }

    // 计算幅度
    static float computeMagnitude(const Complex &sample)
    {
        return std::abs(sample);
    }

    // 计算相位
    static float computePhase(const Complex &sample)
    {
        return std::arg(sample);
    }

    // 汉宁窗
    static std::vector<float> hanningWindow(int size)
    {
        std::vector<float> window(size);
        for (int i = 0; i < size; ++i)
        {
            window[i] = 0.5f * (1.0f - std::cos(2.0f * M_PI * i / (size - 1)));
        }
        return window;
    }

    // 应用窗函数
    static void applyWindow(std::vector<Complex> &data, const std::vector<float> &window)
    {
        for (size_t i = 0; i < data.size(); ++i)
        {
            data[i] *= window[i];
        }
    }

    // 快速傅里叶变换 (使用FFTW库)
    static void fft(std::vector<Complex> &data)
    {
        int n = data.size();
        fftwf_complex *in = reinterpret_cast<fftwf_complex *>(data.data());
        fftwf_complex *out = reinterpret_cast<fftwf_complex *>(data.data());

        fftwf_plan plan = fftwf_plan_dft_1d(n, in, out, FFTW_FORWARD, FFTW_ESTIMATE);
        fftwf_execute(plan);
        fftwf_destroy_plan(plan);
    }

    // 快速傅里叶逆变换
    static void ifft(std::vector<Complex> &data)
    {
        int n = data.size();
        fftwf_complex *in = reinterpret_cast<fftwf_complex *>(data.data());
        fftwf_complex *out = reinterpret_cast<fftwf_complex *>(data.data());

        fftwf_plan plan = fftwf_plan_dft_1d(n, in, out, FFTW_BACKWARD, FFTW_ESTIMATE);
        fftwf_execute(plan);
        fftwf_destroy_plan(plan);

        // 归一化
        for (int i = 0; i < n; ++i)
        {
            data[i] /= static_cast<float>(n);
        }
    }

    // 二维FFT
    static void fft2d(MatrixXcf &data)
    {
        int rows = data.rows();
        int cols = data.cols();

        // 对每行进行FFT
        for (int i = 0; i < rows; ++i)
        {
            std::vector<Complex> row_vec(cols);
            for (int j = 0; j < cols; ++j)
            {
                row_vec[j] = data(i, j);
            }
            fft(row_vec);
            for (int j = 0; j < cols; ++j)
            {
                data(i, j) = row_vec[j];
            }
        }

        // 对每列进行FFT
        for (int j = 0; j < cols; ++j)
        {
            std::vector<Complex> col_vec(rows);
            for (int i = 0; i < rows; ++i)
            {
                col_vec[i] = data(i, j);
            }
            fft(col_vec);
            for (int i = 0; i < rows; ++i)
            {
                data(i, j) = col_vec[i];
            }
        }
    }
};

typedef enum
{
    CA_CFAR, // 单元平均CFAR
    GO_CFAR, // 最大选择CFAR
    SO_CFAR, // 最小选择CFAR
    OS_CFAR  // 有序统计CFAR
} CFARType;

// CFAR检测器
class CFARDetector
{
private:
    int guard_cells;
    int training_cells;
    float threshold_factor;
    CFARType cfar_type;

public:
    CFARDetector(int guard, int train, float factor, CFARType type = CA_CFAR)
        : guard_cells(guard), training_cells(train), threshold_factor(factor), cfar_type(type) {}

    // 检测函数
    std::vector<int> detect(const std::vector<float> &power_spectrum)
    {
        std::vector<int> detections;
        int n = power_spectrum.size();

        for (int i = 0; i < n; ++i)
        {
            float noise_power = 0.0f;

            switch (cfar_type)
            {
            case CA_CFAR:
                noise_power = computeCACFARNoise(power_spectrum, i);
                break;
            case GO_CFAR:
                noise_power = computeGOCFARNoise(power_spectrum, i);
                break;
            case SO_CFAR:
                noise_power = computeSOCFARNoise(power_spectrum, i);
                break;
            case OS_CFAR:
                noise_power = computeOSCFARNoise(power_spectrum, i);
                break;
            }

            // 应用阈值
            if (noise_power > 0 && power_spectrum[i] > threshold_factor * noise_power)
            {
                detections.push_back(i);
            }
        }

        return detections;
    }

private:
    // 单元平均CFAR
    float computeCACFARNoise(const std::vector<float> &power_spectrum, int test_cell)
    {
        int left_start = std::max(0, test_cell - guard_cells - training_cells);
        int left_end = std::max(0, test_cell - guard_cells);

        int right_start = std::min(static_cast<int>(power_spectrum.size()), test_cell + guard_cells + 1);
        int right_end = std::min(static_cast<int>(power_spectrum.size()), test_cell + guard_cells + training_cells + 1);

        float noise_power = 0.0f;
        int count = 0;

        for (int j = left_start; j < left_end; ++j)
        {
            noise_power += power_spectrum[j];
            count++;
        }

        for (int j = right_start; j < right_end; ++j)
        {
            noise_power += power_spectrum[j];
            count++;
        }

        return count > 0 ? noise_power / count : 0.0f;
    }

    // 最大选择CFAR
    float computeGOCFARNoise(const std::vector<float> &power_spectrum, int test_cell)
    {
        int left_start = std::max(0, test_cell - guard_cells - training_cells);
        int left_end = std::max(0, test_cell - guard_cells);

        int right_start = std::min(static_cast<int>(power_spectrum.size()), test_cell + guard_cells + 1);
        int right_end = std::min(static_cast<int>(power_spectrum.size()), test_cell + guard_cells + training_cells + 1);

        float left_noise = 0.0f;
        float right_noise = 0.0f;
        int left_count = 0;
        int right_count = 0;

        for (int j = left_start; j < left_end; ++j)
        {
            left_noise += power_spectrum[j];
            left_count++;
        }

        for (int j = right_start; j < right_end; ++j)
        {
            right_noise += power_spectrum[j];
            right_count++;
        }

        left_noise = left_count > 0 ? left_noise / left_count : 0.0f;
        right_noise = right_count > 0 ? right_noise / right_count : 0.0f;

        return std::max(left_noise, right_noise);
    }

    // 最小选择CFAR
    float computeSOCFARNoise(const std::vector<float> &power_spectrum, int test_cell)
    {
        int left_start = std::max(0, test_cell - guard_cells - training_cells);
        int left_end = std::max(0, test_cell - guard_cells);

        int right_start = std::min(static_cast<int>(power_spectrum.size()), test_cell + guard_cells + 1);
        int right_end = std::min(static_cast<int>(power_spectrum.size()), test_cell + guard_cells + training_cells + 1);

        float left_noise = 0.0f;
        float right_noise = 0.0f;
        int left_count = 0;
        int right_count = 0;

        for (int j = left_start; j < left_end; ++j)
        {
            left_noise += power_spectrum[j];
            left_count++;
        }

        for (int j = right_start; j < right_end; ++j)
        {
            right_noise += power_spectrum[j];
            right_count++;
        }

        left_noise = left_count > 0 ? left_noise / left_count : 0.0f;
        right_noise = right_count > 0 ? right_noise / right_count : 0.0f;

        return std::min(left_noise, right_noise);
    }

    // 有序统计CFAR
    float computeOSCFARNoise(const std::vector<float> &power_spectrum, int test_cell)
    {
        int start = std::max(0, test_cell - guard_cells - training_cells);
        int end = std::min(static_cast<int>(power_spectrum.size()), test_cell + guard_cells + training_cells + 1);

        std::vector<float> window;
        for (int j = start; j < end; ++j)
        {
            if (j < test_cell - guard_cells || j > test_cell + guard_cells)
            {
                window.push_back(power_spectrum[j]);
            }
        }

        if (window.empty())
            return 0.0f;

        // 排序并选择第k个值
        int k = static_cast<int>(window.size() * 0.75f); // 选择75%位置的值
        std::nth_element(window.begin(), window.begin() + k, window.end());

        return window[k];
    }
};

// MTI/MTD处理
class MTI_MTDProcessor
{
private:
    RadarParameters params;
    std::vector<Pulse> delay_line;
    int delay_line_index;

public:
    MTI_MTDProcessor(const RadarParameters &radar_params)
        : params(radar_params), delay_line_index(0)
    {
        // 初始化延迟线
        delay_line.resize(3, Pulse(params.num_samples, Complex(0, 0)));
    }

    // 三脉冲对消器 (MTI)
    std::vector<Pulse> processMTI(const std::vector<Pulse> &pulses)
    {
        std::vector<Pulse> result;
        result.reserve(pulses.size());

        for (const auto &pulse : pulses)
        {
            // 更新延迟线
            delay_line[delay_line_index] = pulse;
            delay_line_index = (delay_line_index + 1) % 3;

            // 三脉冲对消: y[n] = x[n] - 2x[n-1] + x[n-2]
            if (delay_line_index == 0)
            { // 延迟线已填满
                Pulse cancelled_pulse(params.num_samples);
                for (int i = 0; i < params.num_samples; ++i)
                {
                    cancelled_pulse[i] = delay_line[0][i] - 2.0f * delay_line[1][i] + delay_line[2][i];
                }
                result.push_back(cancelled_pulse);
            }
        }

        return result;
    }

    // 脉冲多普勒处理 (MTD)
    PulseDopplerMatrix processMTD(const std::vector<Pulse> &pulses)
    {
        int num_pulses = pulses.size();
        PulseDopplerMatrix doppler_matrix(params.num_samples,
                                          std::vector<Complex>(num_pulses, Complex(0, 0)));

        // 创建窗函数
        auto window = SignalProcessor::hanningWindow(num_pulses);

        // 对每个距离门进行多普勒处理
        for (int range_bin = 0; range_bin < params.num_samples; ++range_bin)
        {
            std::vector<Complex> doppler_profile(num_pulses);

            // 提取同一距离门的所有脉冲
            for (int pulse_idx = 0; pulse_idx < num_pulses; ++pulse_idx)
            {
                doppler_profile[pulse_idx] = pulses[pulse_idx][range_bin];
            }

            // 应用窗函数
            SignalProcessor::applyWindow(doppler_profile, window);

            // 应用FFT进行多普勒分析
            SignalProcessor::fft(doppler_profile);

            // 存储结果
            for (int doppler_bin = 0; doppler_bin < num_pulses; ++doppler_bin)
            {
                doppler_matrix[range_bin][doppler_bin] = doppler_profile[doppler_bin];
            }
        }

        return doppler_matrix;
    }

    // 二维FFT处理 (距离-多普勒)
    PulseDopplerMatrix processRangeDoppler(const std::vector<Pulse> &pulses)
    {
        int num_pulses = pulses.size();

        // 转换为Eigen矩阵以便处理
        MatrixXcf rd_matrix(num_pulses, params.num_samples);

        for (int i = 0; i < num_pulses; ++i)
        {
            for (int j = 0; j < params.num_samples; ++j)
            {
                rd_matrix(i, j) = pulses[i][j];
            }
        }

        // 应用二维FFT
        SignalProcessor::fft2d(rd_matrix);

        // 转换回结果格式
        PulseDopplerMatrix result(params.num_samples,
                                  std::vector<Complex>(num_pulses, Complex(0, 0)));

        for (int i = 0; i < num_pulses; ++i)
        {
            for (int j = 0; j < params.num_samples; ++j)
            {
                result[j][i] = rd_matrix(i, j);
            }
        }

        return result;
    }
};

// 测角处理器 (使用DBF技术)
class AngleEstimator
{
private:
    RadarParameters params;
    std::vector<float> antenna_positions;

public:
    AngleEstimator(const RadarParameters &radar_params)
        : params(radar_params)
    {
        // 计算天线位置 (线性阵列)
        for (int i = 0; i < params.num_antennas; ++i)
        {
            antenna_positions.push_back(i * params.antenna_spacing);
        }
    }

    // 波束形成 - 计算空间谱
    std::vector<float> beamforming(const std::vector<Complex> &antenna_data,
                                   const std::vector<float> &angles)
    {
        std::vector<float> spatial_spectrum(angles.size(), 0.0f);

        for (size_t i = 0; i < angles.size(); ++i)
        {
            float angle_rad = angles[i] * M_PI / 180.0f;
            Complex sum(0, 0);

            for (int j = 0; j < params.num_antennas; ++j)
            {
                // 计算波程差引起的相位差
                float phase_shift = 2.0f * M_PI * antenna_positions[j] *
                                    std::sin(angle_rad) / params.light_speed *
                                    params.carrier_frequency;

                // 应用相位补偿
                sum += antenna_data[j] * std::exp(Complex(0, -phase_shift));
            }

            spatial_spectrum[i] = std::norm(sum);
        }

        return spatial_spectrum;
    }

    // MUSIC算法 - 高分辨率测角
    std::vector<float> musicAlgorithm(const std::vector<Complex> &antenna_data,
                                      const std::vector<float> &angles, int signal_num)
    {
        int snapshot_num = 1; // 假设单次快拍
        MatrixXcf data_matrix(params.num_antennas, snapshot_num);

        // 填充数据矩阵
        for (int i = 0; i < params.num_antennas; ++i)
        {
            data_matrix(i, 0) = antenna_data[i];
        }

        // 计算协方差矩阵
        MatrixXcf covariance = data_matrix * data_matrix.adjoint() / snapshot_num;

        // 特征值分解
        Eigen::SelfAdjointEigenSolver<MatrixXcf> eigensolver(covariance);
        if (eigensolver.info() != Eigen::Success)
        {
            return std::vector<float>(angles.size(), 0.0f);
        }

        // 获取特征值和特征向量
        VectorXcf eigenvalues = eigensolver.eigenvalues();
        MatrixXcf eigenvectors = eigensolver.eigenvectors();

        // 排序特征值 (降序)
        std::vector<int> indices(params.num_antennas);
        for (int i = 0; i < params.num_antennas; ++i)
            indices[i] = i;

        std::sort(indices.begin(), indices.end(),
                  [&eigenvalues](int a, int b)
                  {
                      return std::real(eigenvalues[a]) > std::real(eigenvalues[b]);
                  });

        // 构建噪声子空间
        MatrixXcf noise_subspace(params.num_antennas, params.num_antennas - signal_num);
        for (int i = 0; i < params.num_antennas - signal_num; ++i)
        {
            noise_subspace.col(i) = eigenvectors.col(indices[signal_num + i]);
        }

        // 计算MUSIC谱
        std::vector<float> music_spectrum(angles.size(), 0.0f);

        for (size_t i = 0; i < angles.size(); ++i)
        {
            float angle_rad = angles[i] * M_PI / 180.0f;
            VectorXcf steering_vector(params.num_antennas);

            // 构建导向矢量
            for (int j = 0; j < params.num_antennas; ++j)
            {
                float phase_shift = 2.0f * M_PI * antenna_positions[j] *
                                    std::sin(angle_rad) / params.light_speed *
                                    params.carrier_frequency;
                steering_vector(j) = std::exp(Complex(0, -phase_shift));
            }

            // 计算MUSIC谱
            Complex denominator = steering_vector.adjoint() * noise_subspace *
                                  noise_subspace.adjoint() * steering_vector;
            music_spectrum[i] = 1.0f / std::abs(denominator);
        }

        return music_spectrum;
    }

    // 估计角度
    float estimateAngle(const std::vector<Complex> &antenna_data,
                        const std::vector<float> &angles)
    {
        auto spectrum = beamforming(antenna_data, angles);

        // 找到峰值位置
        auto max_it = std::max_element(spectrum.begin(), spectrum.end());
        int peak_idx = std::distance(spectrum.begin(), max_it);

        return angles[peak_idx];
    }
};

// 点迹凝聚器
class PlotAggregator
{
private:
    float range_threshold;
    float doppler_threshold;
    float angle_threshold;
    float power_threshold;

public:
    PlotAggregator(float range_thresh, float doppler_thresh,
                   float angle_thresh, float power_thresh)
        : range_threshold(range_thresh), doppler_threshold(doppler_thresh),
          angle_threshold(angle_thresh), power_threshold(power_thresh) {}

    // 点迹凝聚算法
    std::vector<TargetInfo> aggregate(const std::vector<DetectionCell> &detections,
                                      const RadarParameters &params)
    {
        if (detections.empty())
        {
            return {};
        }

        std::vector<TargetInfo> targets;
        std::vector<bool> processed(detections.size(), false);

        for (size_t i = 0; i < detections.size(); ++i)
        {
            if (processed[i])
                continue;

            // 创建新目标
            std::vector<size_t> cluster = {i};
            processed[i] = true;

            // 寻找邻近点迹
            for (size_t j = i + 1; j < detections.size(); ++j)
            {
                if (processed[j])
                    continue;

                if (isNeighbor(detections[i], detections[j], params))
                {
                    cluster.push_back(j);
                    processed[j] = true;
                }
            }

            // 计算目标参数
            targets.push_back(computeTargetParameters(cluster, detections, params));
        }

        return targets;
    }

private:
    // 判断两个点迹是否相邻
    bool isNeighbor(const DetectionCell &a, const DetectionCell &b,
                    const RadarParameters &params)
    {
        float range_diff = std::abs(calculateRange(a.range_bin, params) -
                                    calculateRange(b.range_bin, params));
        float velocity_diff = std::abs(calculateVelocity(a.doppler_bin, params) -
                                       calculateVelocity(b.doppler_bin, params));
        float angle_diff = std::abs(calculateAngle(a.angle_bin, params) -
                                    calculateAngle(b.angle_bin, params));

        return range_diff <= range_threshold &&
               velocity_diff <= doppler_threshold &&
               angle_diff <= angle_threshold;
    }

    // 计算目标参数
    TargetInfo computeTargetParameters(const std::vector<size_t> &cluster,
                                       const std::vector<DetectionCell> &detections,
                                       const RadarParameters &params)
    {
        float range_sum = 0.0f;
        float velocity_sum = 0.0f;
        float angle_sum = 0.0f;
        float power_sum = 0.0f;
        float max_power = -std::numeric_limits<float>::max();

        for (size_t idx : cluster)
        {
            const auto &det = detections[idx];
            float range = calculateRange(det.range_bin, params);
            float velocity = calculateVelocity(det.doppler_bin, params);
            float angle = calculateAngle(det.angle_bin, params);

            range_sum += range;
            velocity_sum += velocity;
            angle_sum += angle;
            power_sum += det.power;
            max_power = std::max(max_power, det.power);
        }

        int count = cluster.size();
        float avg_range = range_sum / count;
        float avg_velocity = velocity_sum / count;
        float avg_angle = angle_sum / count;
        float avg_power = power_sum / count;

        // 计算SNR (简化计算)
        float noise_floor = avg_power - 10.0f; // 假设噪声比平均功率低10dB
        float snr = max_power - noise_floor;

        return TargetInfo(avg_range, avg_velocity, avg_angle, avg_power, snr, count);
    }

    // 计算距离
    float calculateRange(int range_bin, const RadarParameters &params)
    {
        return range_bin * params.getRangeResolution();
    }

    // 计算速度
    float calculateVelocity(int doppler_bin, const RadarParameters &params)
    {
        float doppler_freq = (doppler_bin - params.num_pulses / 2) *
                             (1.0f / (params.num_pulses * params.pulse_rep_interval));
        return doppler_freq * params.light_speed / (2.0f * params.carrier_frequency);
    }

    // 计算角度
    float calculateAngle(int angle_bin, const RadarParameters &params)
    {
        return static_cast<float>(angle_bin) * 180.0f / params.num_antennas - 90.0f;
    }
};

// 主雷达处理器
class RadarProcessor
{
private:
    RadarParameters params;
    std::unique_ptr<CFARDetector> cfar_detector;
    std::unique_ptr<MTI_MTDProcessor> mti_mtd_processor;
    std::unique_ptr<AngleEstimator> angle_estimator;
    std::unique_ptr<PlotAggregator> plot_aggregator;

public:
    RadarProcessor(const RadarParameters &radar_params)
        : params(radar_params)
    {
        // 初始化各处理器
        cfar_detector = std::make_unique<CFARDetector>(2, 8, 1.5f, CA_CFAR); // CFARDetector::
        mti_mtd_processor = std::make_unique<MTI_MTDProcessor>(params);
        angle_estimator = std::make_unique<AngleEstimator>(params);
        plot_aggregator = std::make_unique<PlotAggregator>(50.0f, 5.0f, 5.0f, 3.0f);
    }

    // 处理雷达数据
    std::vector<TargetInfo> process(const std::vector<std::vector<Pulse>> &antenna_data)
    {
        std::vector<DetectionCell> all_detections;

        // 对每个天线通道进行处理
        for (int ant_idx = 0; ant_idx < params.num_antennas; ++ant_idx)
        {
            // 1. MTI/MTD处理
            auto mti_data = mti_mtd_processor->processMTI(antenna_data[ant_idx]);
            auto mtd_matrix = mti_mtd_processor->processMTD(mti_data);

            // 2. 计算功率谱
            std::vector<std::vector<float>> power_spectrum(
                params.num_samples,
                std::vector<float>(params.num_pulses - 1, 0.0f));

            for (int i = 0; i < params.num_samples; ++i)
            {
                for (int j = 0; j < params.num_pulses - 1; ++j)
                {
                    power_spectrum[i][j] = SignalProcessor::computePowerDB(mtd_matrix[i][j]);
                }
            }

            // 3. CFAR检测
            for (int i = 0; i < params.num_samples; ++i)
            {
                auto detections = cfar_detector->detect(power_spectrum[i]);
                for (int doppler_bin : detections)
                {
                    DetectionCell cell;
                    cell.range_bin = i;
                    cell.doppler_bin = doppler_bin;
                    cell.angle_bin = ant_idx; // 暂时存储天线索引
                    cell.power = power_spectrum[i][doppler_bin];
                    all_detections.push_back(cell);
                }
            }
        }

        // 4. 角度估计
        estimateAngles(all_detections, antenna_data);

        // 5. 点迹凝聚
        return plot_aggregator->aggregate(all_detections, params);
    }

private:
    // 角度估计
    void estimateAngles(std::vector<DetectionCell> &detections,
                        const std::vector<std::vector<Pulse>> &antenna_data)
    {
        // 生成角度搜索范围
        std::vector<float> angles(181);
        for (int i = 0; i <= 180; ++i)
        {
            angles[i] = i - 90.0f; // -90° 到 +90°
        }

        for (auto &det : detections)
        {
            // 提取同一距离-多普勒单元的所有天线数据
            std::vector<Complex> ant_samples(params.num_antennas);
            for (int ant_idx = 0; ant_idx < params.num_antennas; ++ant_idx)
            {
                // 注意：这里需要根据MTD处理后的数据索引来获取
                // 简化处理，实际应根据具体数据结构调整
                ant_samples[ant_idx] = antenna_data[ant_idx][det.doppler_bin][det.range_bin];
            }

            // 估计角度
            det.angle_bin = static_cast<int>(angle_estimator->estimateAngle(ant_samples, angles) + 90.0f);
        }
    }
};

// 数据模拟器
class DataSimulator
{
private:
    RadarParameters params;
    std::default_random_engine generator;
    std::normal_distribution<float> noise_dist;

public:
    DataSimulator(const RadarParameters &radar_params, float noise_stddev = 0.1f)
        : params(radar_params), noise_dist(0.0f, noise_stddev) {}

    // 生成模拟雷达数据
    std::vector<std::vector<Pulse>> generateData(const std::vector<TargetInfo> &targets)
    {
        std::vector<std::vector<Pulse>> antenna_data(
            params.num_antennas,
            std::vector<Pulse>(params.num_pulses,
                               Pulse(params.num_samples, Complex(0, 0))));

        // 添加噪声
        for (int ant = 0; ant < params.num_antennas; ++ant)
        {
            for (int pulse = 0; pulse < params.num_pulses; ++pulse)
            {
                for (int sample = 0; sample < params.num_samples; ++sample)
                {
                    float real = noise_dist(generator);
                    float imag = noise_dist(generator);
                    antenna_data[ant][pulse][sample] = Complex(real, imag);
                }
            }
        }

        // 添加目标信号
        for (const auto &target : targets)
        {
            addTarget(antenna_data, target);
        }

        return antenna_data;
    }

private:
    void addTarget(std::vector<std::vector<Pulse>> &data, const TargetInfo &target)
    {
        // 计算目标参数对应的单元索引
        int range_bin = static_cast<int>(target.range / params.getRangeResolution());
        float doppler_freq = 2.0f * target.velocity * params.carrier_frequency / params.light_speed;
        int doppler_bin = static_cast<int>(doppler_freq * params.num_pulses * params.pulse_rep_interval) +
                          params.num_pulses / 2;

        float angle_rad = target.angle * M_PI / 180.0f;

        // 添加目标回波
        for (int ant = 0; ant < params.num_antennas; ++ant)
        {
            // 计算天线间的相位差
            float ant_phase = 2.0f * M_PI * ant * params.antenna_spacing *
                              std::sin(angle_rad) / params.light_speed *
                              params.carrier_frequency;

            for (int pulse = 0; pulse < params.num_pulses; ++pulse) // -----------------------------------------------------------------------params.num_pulses
            {
                // 计算多普勒相位
                float doppler_phase = 2.0f * M_PI * doppler_freq * pulse * params.pulse_rep_interval;

                // 计算信号幅度 (线性)
                float amplitude = std::sqrt(std::pow(10.0f, target.power / 10.0f));

                // 生成目标信号
                Complex signal = amplitude * std::exp(Complex(0, doppler_phase + ant_phase));

                // 添加到数据中
                if (range_bin >= 0 && range_bin < params.num_samples)
                {
                    data[ant][pulse][range_bin] += signal;
                }
            }
        }
    }
};

// 示例用法
int main()
{
    // 设置雷达参数
    RadarParameters params;
    params.carrier_frequency = 10e9f;  // 10 GHz
    params.pulse_width = 1e-6f;        // 1 μs
    params.pulse_rep_interval = 1e-3f; // 1 ms
    params.num_pulses = 16;            // 64;
    params.num_samples = 1024;
    params.sampling_rate = 10e6f; // 10 MHz
    params.bandwidth = 5e6f;      // 5 MHz
    params.num_antennas = 8;
    params.antenna_spacing = 0.015f; // 15 mm
    params.light_speed = 3e8f;       // 光速

    std::cout << "雷达参数:" << std::endl;
    std::cout << "距离分辨率: " << params.getRangeResolution() << " m" << std::endl;
    std::cout << "最大不模糊距离: " << params.getMaxUnambiguousRange() << " m" << std::endl;
    std::cout << "速度分辨率: " << params.getVelocityResolution() << " m/s" << std::endl;
    std::cout << "最大不模糊速度: " << params.getMaxUnambiguousVelocity() << " m/s" << std::endl;
    std::cout << std::endl;

    // 创建雷达处理器
    RadarProcessor radar_processor(params);

    // 创建数据模拟器
    DataSimulator data_simulator(params);

    // 定义模拟目标
    std::vector<TargetInfo> simulated_targets = {
        {5000.0f, 100.0f, 30.0f, 60.0f, 0.0f, 0}, // 5km, 100m/s, 30°, 60dB
                                                  // {8000.0f, -50.0f, -10.0f, 55.0f, 0.0f, 0}, // 8km, -50m/s, -10°, 55dB
        //{3000.0f, 0.0f, 5.0f, 50.0f, 0.0f, 0}      // 3km, 0m/s, 5°, 50dB
    };

    // 生成模拟数据
    auto simulated_data = data_simulator.generateData(simulated_targets);

    // 处理数据
    auto detected_targets = radar_processor.process(simulated_data);

    // 输出结果
    std::cout << "检测到的目标:" << std::endl;
    int i = 0;
    for (const auto &target : detected_targets)
    {
        std::cout << "距离: " << target.range << " m, "
                  << "速度: " << target.velocity << " m/s, "
                  << "角度: " << target.angle << " deg, "
                  << "功率: " << target.power << " dB, "
                  << "SNR: " << target.snr << " dB, "
                  << "点迹数: " << target.num_detections << std::endl;

        i++;
    }
    std::cout << "total point:" << i << std::endl;
    return 0;
}